## Function to install and load packages
loadPkg=function(toLoad){
  for(lib in toLoad){
    if(! lib %in% installed.packages()[,1])
    { install.packages(lib, repos='http://cran.rstudio.com/') }
    suppressMessages( library(lib, character.only=TRUE) ) }
}


# hurdle model SE
clx2 <- 
  function(fm, dfcw,  cluster){
    # R-codes (www.r-project.org) for computing
    # clustered-standard errors. Mahmood Arai, Jan 26, 2008.
    
    # The arguments of the function are:
    # fitted model, cluster1
    # You need to install libraries `sandwich' and `lmtest'
    
    # reweighting the var-cov matrix for the within model
    require(sandwich);require(lmtest)
    if(class(fm)=="zeroinfl"|class(fm)=="hurdle") {
      M <- length(unique(cluster))   
      N <- length(cluster)           
      K <- dim(fm$vcov)[1]        #here is the rank from the zero inflated fits             
      dfc <- (M/(M-1))*((N-1)/(N-K))  
      uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
      vcovCL <- dfc[1]*sandwich(fm, meat=crossprod(uj)/N)*dfcw #fix a length problem in dfc
      list(summary=coeftest(fm, vcovCL))}
    else if(class(fm)!="zeroinfl"){
      M <- length(unique(cluster))
      N <- length(cluster)
      K <- fm$rank
      dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
      uj <- apply(estfun(fm), 2, function(x) tapply(x, cluster, sum));
      rcse.cov <- dfc * sandwich(fm, meat = crossprod(uj)/N)
      rcse.se <- coeftest(fm, rcse.cov)
      return(list( rcse.se))}
    
  }





# clustered SE
# function that does clustered SEs
vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}


source('/Users/JuamnTellez/Dropbox/github/stargazer/R/stargazer-internal.R')
stargazer <-
  function(..., type = "latex", title="", style="default", summary=NULL, out=NULL, out.header=FALSE,
           column.labels=NULL, column.separate = NULL, covariate.labels=NULL, dep.var.caption=NULL, 
           dep.var.labels=NULL, dep.var.labels.include=TRUE, align=FALSE, coef=NULL, se=NULL, t=NULL, 
           p=NULL,  t.auto=TRUE, p.auto=TRUE, ci=FALSE, ci.custom=NULL, ci.level=0.95, ci.separator=NULL, 
           add.lines=NULL, apply.coef=NULL, apply.se=NULL, apply.t=NULL, apply.p=NULL, apply.ci=NULL, 
           colnames = NULL, column.sep.width = "5pt", 
           decimal.mark=NULL, df=TRUE, digit.separate=NULL, digit.separator=NULL, digits=NULL, digits.extra=NULL, 
           flip=FALSE,
           float=TRUE, float.env="table", font.size=NULL, header=TRUE, initial.zero=NULL, intercept.bottom=TRUE, 
           intercept.top=FALSE, keep=NULL, keep.stat=NULL, label="", model.names=NULL, model.numbers=NULL, 
           multicolumn=TRUE, no.space=NULL, notes=NULL, notes.align=NULL, notes.append=TRUE, notes.label=NULL, 
           object.names=FALSE,
           omit=NULL, omit.labels=NULL, omit.stat=NULL, omit.summary.stat=NULL, omit.table.layout=NULL,
           omit.yes.no=c("Yes","No"), order=NULL, ord.intercepts=FALSE, adjust.box = TRUE,
           perl=FALSE, report=NULL, rownames = NULL, 
           rq.se = "nid", selection.equation=FALSE, single.row=FALSE, star.char=NULL, 
           star.cutoffs=NULL, suppress.errors=FALSE, table.layout=NULL, table.placement = "!htbp", zero.component=FALSE, 
           summary.logical=TRUE, summary.stat=NULL, nobs=TRUE, mean.sd=TRUE, min.max=TRUE, median=FALSE, 
           iqr=FALSE) {
    
    save.warn.option <- getOption("warn") 
    options(warn=-1)
    return(.stargazer.wrap(..., type=type, title=title, style=style, summary=summary, out=out, out.header=out.header,
                           column.labels=column.labels, column.separate = column.separate, adjust.box = adjust.box,
                           covariate.labels=covariate.labels, dep.var.caption = dep.var.caption,
                           dep.var.labels=dep.var.labels, dep.var.labels.include=dep.var.labels.include,
                           align=align, coef=coef, se=se, t=t, p=p, t.auto=t.auto, p.auto=p.auto, 
                           ci=ci, ci.custom=ci.custom, ci.level=ci.level, ci.separator = ci.separator,
                           add.lines=add.lines, apply.coef=apply.coef, apply.se=apply.se, apply.t=apply.t, 
                           apply.p=apply.p, apply.ci=apply.ci, colnames=colnames,
                           column.sep.width=column.sep.width, decimal.mark=decimal.mark, df=df,
                           digit.separate=digit.separate, digit.separator=digit.separator, 
                           digits=digits, digits.extra=digits.extra, 
                           flip=flip,
                           float=float, float.env=float.env,
                           font.size=font.size, header=header, 
                           initial.zero=initial.zero, 
                           intercept.bottom=intercept.bottom, intercept.top=intercept.top, 
                           keep = keep, keep.stat = keep.stat,
                           label = label,
                           model.names=model.names, model.numbers=model.numbers,
                           multicolumn = multicolumn,
                           no.space=no.space, notes=notes, notes.align=notes.align, 
                           notes.append=notes.append, notes.label=notes.label, object.names=object.names,
                           omit=omit, omit.labels=omit.labels, omit.stat=omit.stat, omit.summary.stat=omit.summary.stat,
                           omit.table.layout=omit.table.layout,omit.yes.no=omit.yes.no, 
                           order=order, ord.intercepts=ord.intercepts, perl=perl, 
                           report=report, rownames=rownames, rq.se=rq.se, selection.equation=selection.equation,
                           single.row=single.row, star.char=star.char, 
                           star.cutoffs=star.cutoffs, suppress.errors=suppress.errors,
                           table.layout=table.layout,
                           table.placement = table.placement, zero.component=zero.component,
                           summary.logical = summary.logical,
                           summary.stat = summary.stat,
                           nobs=nobs, mean.sd=mean.sd, 
                           min.max=min.max, median=median, iqr=iqr, warn=save.warn.option))
  }



## ggcoef
ggCoef = function(modelframe, subset = F)
{
  # Specify the width of your confidence intervals
  interval1 = -qnorm((1-0.9)/2)  # 90% multiplier
  interval2 = -qnorm((1-0.95)/2)  # 95% multiplier
  
  if(subset == F) {
    zp1 = ggplot(modelframe)
    zp1 = zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
    zp1 = zp1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                   ymax = Coefficient + SE*interval1),
                               lwd = 1, position = position_dodge(width = 1/2))
    zp1 = zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                    ymax = Coefficient + SE*interval2),
                                lwd = 1/2, position = position_dodge(width = 1/2),
                                shape = 21, fill = "WHITE")
    zp1 = zp1 + coord_flip() + labs(x = "")
    zp1 = zp1 + theme_bw(base_family = 'Helvetica Neue')
    zp1 = zp1 + theme(axis.text.x = element_text(family = 'Helvetica Neue Light'), 
                      axis.text.y = element_text(family = 'Helvetica Neue Light'), 
                      axis.ticks = element_blank(), 
                      panel.border = element_blank(), 
                      legend.background=element_blank())
    return(zp1)
  } else {
    zp1 = ggplot(modelframe, aes(colour = Subset, shape = Subset))
    zp1 = zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
    zp1 = zp1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                   ymax = Coefficient + SE*interval1),
                               lwd = 1, position = position_dodge(width = 1/2))
    zp1 = zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                    ymax = Coefficient + SE*interval2),
                                lwd = 1/2, position = position_dodge(width = 1/2), fill = "WHITE")
    zp1 = zp1 + coord_flip() + labs(x = "", colour = 'Sample subset')
    zp1 = zp1 + theme_bw(base_family = 'Helvetica Neue')
    zp1 = zp1 + theme(axis.text.x = element_text(family = 'Helvetica Neue Light'), 
                      axis.text.y = element_text(family = 'Helvetica Neue Light'), 
                      axis.ticks = element_blank(), 
                      panel.border = element_blank(), 
                      legend.background=element_blank())
    return(zp1)
  }
}

